home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / mathpack / cffti1.f < prev    next >
Encoding:
Text File  |  1989-08-14  |  1.7 KB  |  73 lines
  1. *deck cffti1
  2.       subroutine cffti1 (n,wa,ifac)
  3. C***BEGIN PROLOGUE  CFFTI1
  4. C***REFER TO CFFTI
  5. C***ROUTINES CALLED  (NONE)
  6. C***END PROLOGUE  CFFTI1
  7.       dimension       wa(1)      ,ifac(1)    ,ntryh(4)
  8.       data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/
  9. C***FIRST EXECUTABLE STATEMENT  CFFTI1
  10.       nl = n
  11.       nf = 0
  12.       j = 0
  13.   101 j = j+1
  14.       if (j-4) 102,102,103
  15.   102 ntry = ntryh(j)
  16.       go to 104
  17.   103 ntry = ntry+2
  18.   104 nq = nl/ntry
  19.       nr = nl-ntry*nq
  20.       if (nr) 101,105,101
  21.   105 nf = nf+1
  22.       ifac(nf+2) = ntry
  23.       nl = nq
  24.       if (ntry .ne. 2) go to 107
  25.       if (nf .eq. 1) go to 107
  26.       do 106 i=2,nf
  27.          ib = nf-i+2
  28.          ifac(ib+2) = ifac(ib+1)
  29.   106 continue
  30.       ifac(3) = 2
  31.   107 if (nl .ne. 1) go to 104
  32.       ifac(1) = n
  33.       ifac(2) = nf
  34.       tpi = 6.28318530717959
  35.       argh = tpi/float(n)
  36.       i = 2
  37.       l1 = 1
  38.       do 110 k1=1,nf
  39.          ip = ifac(k1+2)
  40.          ld = 0
  41.          l2 = l1*ip
  42.          ido = n/l2
  43.          idot = ido+ido+2
  44.          ipm = ip-1
  45.          wldc = 1.
  46.          wlds = 0.
  47.          arg = float(l1)*argh
  48.          dl1c = cos(arg)
  49.          dl1s = sin(arg)
  50.          do 109 j=1,ipm
  51.             i1 = i
  52.             wa(i-1) = 1.
  53.             wa(i) = 0.
  54.             ld = ld+l1
  55.             wldch = wldc
  56.             wldc = dl1c*wldc-dl1s*wlds
  57.             wlds = dl1s*wldch+dl1c*wlds
  58.             dc = wldc
  59.             ds = wlds
  60.             do 108 ii=4,idot,2
  61.                i = i+2
  62.                wa(i-1) = dc*wa(i-3)-ds*wa(i-2)
  63.                wa(i) = ds*wa(i-3)+dc*wa(i-2)
  64.   108       continue
  65.             if (ip .le. 5) go to 109
  66.             wa(i1-1) = wa(i-1)
  67.             wa(i1) = wa(i)
  68.   109    continue
  69.          l1 = l2
  70.   110 continue
  71.       return
  72.       end
  73.